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Abstract. We study the pair contact process with diffusion (PCPD) using Monte 
Carlo simulations, and concentrate on the decay of the particle density p with time, 
near its critical point, which is assumed to follow p(t) « ct~ s + c 2 t~ S2 + . . .. This model 
is known for its slow convergence to the asymptotic critical behavior; we therefore pay 
particular attention to finite-time corrections. We find that at the critical point, the 
ratio of p and the pair density p p converges to a constant, indicating that both densities 
decay with the same powerlaw. We show that under the assumption 62 ~ 2<5, two of the 
critical exponents of the PCPD model are 6 — 0.165(10) and f3 = 0.31(4), consistent 
with those of the directed percolation (DP) model. 
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1. Introduction 

The pair contact process with diffusion (PCPD) is a one-dimensional model of fermionic 
particles on a lattice. In this context, fermionic means that a site cannot be occupied 
by more than one particle. When two particles are adjacent, they can interact in two 
ways: they can annihilate each other, or alternatively, they can create another particle 
on an adjacent lattice site. Particles can also diffuse by hopping from one site to the 
next. The reactions of the PCPD model and their rates are given by: 

I " A Ai -ch with «*» 

[ OAA AAA 2 

AA 00 with rate p (1 - d) ^ ' 

AO -H- OA with rate d 

Given a value for the diffusion coefficient d, we can discern three different regimes 
depending on the annihilation rate p. If p is very large, then annihilation dominates 
the process, and the particles on the lattice will die out quickly. This is the inactive 
phase. On the other hand, if p is very small, the particle creation reaction will ensure 
that (with extremely high probability) the system will maintain a high density. This 
is called the active phase. Well into the active or the inactive phase, long-ranged 
interactions are absent, and both these regimes can therefore be described well by 
mean-field theory. However, if the boundary between the active and inactive phase 
is approached from either side, length and time scales diverge in a power-law fashion. 
In analogy with equilibrium phase transitions, one expects that the system then exhibits 
critical behaviour in the transition between these regimes, at the critical value p c . 
From equilibrium statistical physics it is well established that phase transitions can 
be classified in universality classes, each of which characterized by a unique set of 
critical exponents. Moreover, these exponents are typically insensitive to small changes 
in the model, such as details of short-ranged interactions. A central question in 
non-equilibrium statistical physics is whether also dynamical phase transitions can be 
classified into universality classes, which also are insensitive to such small changes. 

In this paper we will concentrate on the scaling relations concerning the particle 
density p, which are given by: 

Pp= Pc ~ t~ S (e = 0) 

Pt-oo ~ (e > 0) ' 1 ] 

where e = \p — p c \ is the distance from criticality, and 5 and /3 are two of the critical 
exponents of the PCPD system. A conjecture by Grassberger [1] and Janssen [2] states 
that all systems with a single order parameter and a single absorbing state will belong 
to the universality class of the Directed Percolation (DP) model. The critical exponents 
of the one-dimensional DP model are known to very high accuracy (S = 0.159464(6) and 
P = 0.276486(8)) [3]. The values for 5 and of the PCPD model have been disputed 
extensively, with estimates of 5 ranging from 0.16(1) to 0.27(4) and estimates of f3 



Critical exponents of the pair contact process with diffusion 



3 



varying from (3 < 0.34 to = 0.58(1) H El El El El El 101 EU H21 EU El CG2 ESI EH EE] • For 

a more extensive overview, see [3]. The main interest in the PCPD model is whether it 
shares its critical exponents with DP, or not; the latter would disprove the Grassberger- 
Jansen conjecture. 

The difficulties with the determination of the exponent 5 of the PCPD model is 
that there are large finite-time corrections [HES], and the effective exponent thus shows 
a drift with simulation time. In this paper we generate high-quality data, exploiting 
the computational power of graphics processing units (GPUs); a description of the 
implementation will be published elsewhere. We then analyze these data in a way that 
suppresses the leading finite-time corrections. 



2. Simulation results 



At criticality, the asymptotic decay of the density is described by a power-law of the 
form p{t) ~ t" s . For the DP model the corrections to this power-law decay very rapidly, 
and therefore <5dp is known with very high accuracy. However, as the history of the 
PCPD model shows, the density decay in the PCPD model does not show such a clean 
power-law in the time range accessible to computer simulations. Consequently, one 
must account for finite-time corrections, in order to obtain an accurate estimate of S. 
Different finite-time corrections have been proposed in the past, including logarithmic 
corrections [H] and power- law corrections [13]. Our data suggest power- law corrections, 
so that the density decay at the critical point is given by 

pit) = Cl t~ s + c 2 r S2 + ..., (3) 

with 6 < 62 < ■ ■ ■■ By differentiating the logarithm of the density vs. the logarithm 
of the time, one can define an effective decay exponent as 

alog(t) ci 

Note that in the limit of infinite time, the effective exponent <5 e fr approaches the true 
asymptotic value 5. Moreover, the direction from which it approaches the asymptotic 
value is governed by the sign of c 2 /ci, and the speed of convergence by the gap 5 2 — 8. 

At this point, we want to make some practical remarks: 

i) In practice with simulation data, the differentiation is carried out numerically, 
and one makes the approximation 5 c g ~ ~ l °s(A 2t )/ p(t)) ^ 

ii) When simulation measurements of 5 e s are presented, it is usually more convenient 
for the eye if the presentation is such that the infinite-time behavior is inside the plot, 
rather than at infinite distance. A common approach to achieve this, is to plot 5 c s as 
a function of 1/t; this is not a good idea! If the gap is less than 1, the approach to the 
vertical axis will eventually become infinitely steep, and it is hard to predict how far 
an infinitely steep curve will shoot up (or down). Ideally, one should therefore plot 8 e s 
as a function of t~( S2 ~ s \ as then the approach to the vertical axis will follow a straight 
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line (and indirectly provide information on 82 — 8). Alternatively, one can plot 8 e s as a 
function of p^/ 6 ~ l ). 

iii) Besides the particle density p, another numerically accessible quantity is the 
pair density p p , defined as the fraction of neighboring sites, both of which are being 
occupied by a particle. For the pair density at the critical point, relations equivalent to 
eqs.([3]) and (jl]) can be defined. Barkema and Carlon [13] provided numerical evidence 
that the ratio of the densities of single particles and pairs tends to a constant; further on 
in this paper we will confirm this. This has as a consequence that the asymptotic density 
decay of singles and pairs is governed by a unique exponent 8. Also, this means that as 
long as the effective exponents S c g and S P)C g do not coincide, finite-time corrections are 
still significant. 

The first results, presented in figure [U are measurements of 8 c g as a function 
of p, for values of p close to the critical point, on either side. The figure shows 
that 5 c ff approaches its asymptotic value from above. Assuming that at the end of 
our simulations, corrections are dominated by 82, this indicates that C2 is positive. 
Furthermore, under this assumption, the maximal value of 5 c g with a p value lying in 
the inactive phase provides an upper bound for the asymptotic exponent: 5 < 0.19 
at p = 0.15247. This bound is less tight then the bound reported by Hinrichsen|16j. 
because it depends strongly on the maximal simulation time of the dataset, which was 
higher in Hinrichsen's simulation, and less on the accuracy of the dataset, which is 
higher in our simulations. We do however confirm the conclusion of Hinrichsen that 
most of the values reported for the exponent 5 in PCPD are ruled out by this upper 
bound. 

Interestingly, the data in figure [T] are consistent with a linear convergence to 
an asymptotic value of 5 = 5 dp- The most important conclusion drawn from this 
observation is that our data do not provide evidence of a violation of the Grassberger- 
Janssen conjecture; and since the product of system size, length of the simulations, and 
the number of these is higher in our simulations than in earlier studies which claim to 
provide evidence for violation of the conjecture, our conclusion is that earlier claims 
of such violation based on numerical measurements of 5 e g are ill-founded. A second 
observation is that the observed linear convergence indicates that 82/ 5—1 ~ 1 and thus 
that 82 ~ 28. Although our simulations are more extensive than those of earlier reports, 
the data is still not accurate enough to unambiguously rule out combinations of (8, 8 2 ) 
which are slightly different from (8dp, 25dp)', in particular there is quite some ambiguity 
in 8 2 . 

We now turn to an accurate estimation of the critical point p c . For this, we want 
to capitalize on our numerical knowledge on the value for the correction exponent 82, 
to suppress finite-time corrections. We do so by defining an adjusted density: 

p*(t) = 2 P (t)-p(2^ s H), (5) 

for some assumed value of 82. Asymptotically, the densities pit) and p*{t) at the 
critical point will coincide, but at long but not infinite times, the corrected density 
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Figure 1. The effective exponent 5 e g as obtained for both the particle density p 
(colored black) and the pair density p p (gray), as a function of the particle density, with 
a diffusion factor d = 0.5. For d = 0.25 and d = 0.75, we found similar results. Starting 
from low to high in the picture, the simulations are performed with p = 0.15246 (1600 
runs), 0.15247 (3300 runs), 0.152475 (7400 runs) and 0.152485 (1000 runs). The lattice 
contained L = 2 18 = 262144 sites. Each simulation reached a time of approximately 
t = 1.6 ■ 10 7 . The curve for p = 0.152475 is closest to the critical value, which we 
estimate to be p c = 0.152473(2). The straight line signifies a possible extrapolation to 
S = S DP . 

suppresses correction terms with exponents Si close to 62 ; with the correct numerical 
choice for <5 2 the first correction term will even be completely removed, and for a 
reasonable approximation of 62, the power-law decay of p*{t) should be much cleaner 
than that of p(t). 

Figure [2] shows the corrected density and corrected pair density as a function of 
time in a log-log plot, with the choice 62 = 2<5dp = 0.319. This gives us the estimate 
5 = 0.165(10). We conclude that the choice of 62 tightly defines our estimates of both 
5 and p c . 

An accurate method to determine the critical point p c is to use a data collapse 
method on the corrected density p*. The data collapse method maps the dataset with 
off-critical p values onto a single curve for p > p c and a second curve for p < p c , if we 
select the correct values for 5, f3 and p c . To achieve this, the (pair) density and the time 
are rescaled as follows: 

t' = e^ s t, p' = e-y . (6) 

Using the diffusion factor d = 0.5, we obtain a very good data collapse, shown 
in figure El We find that the critical value is p c = 0.152473(2). Our estimate of the 
exponent /3 is less precise then our estimate of 6: (3 = 0.31(4), which is consistent with 
the value known for DP: = 0.2765. The fact that the curve for p = 0.15240 shows a 
statistically significant increase, after reaching a (pseudo-)equilibrium, indicates that the 
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Figure 2. Corrected particle density p* (black) and corrected pair density p* (gray) 
as a function of time, for d = 0.25, p = 0.125142 (dashed); d = 0.5, p = 0.152475 
(solid); and d = 0.75, p = 0.191790 (dotted-dashed). In the correction procedure, we 
used 62 = 0.319 and a lattice size of L — 2 18 = 262144. The thick black line shows a 
possible line with 5 = <5dp- 

corrections due to finite time and due to off-critical values for p are not simply additive: 
our procedure to suppress the leading finite-time corrections at p c is not working equally 
well for the off-critical curves. The 'best' value we find for /3 is therefore still experiencing 
significant corrections to scaling. 

The question that remains is how these parameters depend on our choice for 5 2 - 
We used the same methods described above to obtain estimates with 62 ranging from 
0.319 to 0.5. The results are presented in Table [TJ Note that for our choice 62 = 0.319, 
the leading exponent 5 shows little dependence on the diffusion coefficient d. 

In Ref. jl] the leading correction term ~ t~ S2 was removed by using a linear 
combination of the density p and the pair density p* . This is a valid strategy as long 
as the dominant correction term (c2t~ 52 ) does not cancel in the ratio p/p p . In that case 
this ratio will asymptotically go to a constant with a power law (at criticality): 

p/p p = h + k 2 t s ^ + ... (7) 

Thus, if the density ratio p/p p is plotted against t 5 ~ 5i , one should find a straight 
line (with finite-time corrections from higher-order correction terms). Since we found 
previously that 82 ~ 25, we have plotted the ratio versus the density p in figure HI The 
density ratios for the three values of d are clearly not approaching the vertical axis under 
a fixed angle in this plot, but seem to arrive horizontally at p = 0. Thus, we find that 
the leading correction disappears in the ratio and the method described in [3] fails. 

We plotted the ratio against t s ~ A and we found that the most straight line has a 
value of A 0.51. The large discrepancy between 82 and A increases our confidence in 
the fact that the correction term ~ t~ &2 vanishes in the ratio p/p p . 
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Figure 3. Data collapse of curves obtained with various values for p at d = 0.5. 
The rescaled corrected density is plotted as a function of rescaled time (see Eqs. ([5])). 
The data are obtained from the same simulations as in Fig. [T] added to these data are 
simulations withp = 0.1524 (500 runs on a smaller lattice with L = 65536), p = 0.15242 
(1000 runs), p = 0.15245 (1000 runs), p = 0.15248 (500 runs) and p = 0.1525 (1000 
runs). In the correction procedure, we used 82 = 0.319. In the rescaling, we used 
5 = 0.159 and /3 = 0.2765 as in the DP model, combined with p c = 0.152473. 



Table 1. Values for the exponents 8 and (3 and critical point p c , obtained by our 
analysis approach, for three values of the diffusion coefficient d and for our estimated 
value 82 — 0.319, as well as for some other values for 82- 
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Figure 4. Left panel: ratio p/p p of the particle and pair densities, as a function 
of the particle density. Right panel: same data, plotted as a function of t s ~ A with 
A = 0.51. The different curves correspond to d = 0.25 (solid line), d = 0.5 (larger 
dashes) and d = 0.75 (smaller dashes). Curves are shifted vertically by some arbitrary, 
d-dependent constant c(d), to let them fit nicely into a single figure; lower curves 
correspond to higher values of d. 

Thus, the ratio between the particle density and the pair density goes to a constant 
even faster than expected. This still means that the exponent 8 is equal for the particle 
and the pair density. 

3. Summary and conclusion 

We have performed extensive simulations of the one-dimensional PCPD model, using a 
highly efficient GPU-based simulation approach. The analysis of the simulation results 
was performed in a way that takes account of finite-time effects. Our main goal was to 
verify the Grassberger-Janssen conjecture, which predicts that the exponents for PCPD 
coincide with those of directed percolation. 

We find that our data are consistent with DP values for the exponents 8 which 
describes the decay of the particle density at the critical point, and /3 which describes 
the asymptotic particle density for simulations close to the critical point, but slightly in 
the active phase. 

Additionally, we find that the the leading correction exponent 82 is numerically 
found to be close to 28, which would suggest corrections to scaling of order p 2 . We also 
find that these leading corrections do not seem to manifest themselves in the ratio of 
the particle and pair densities. 

Under the assumption that the Grassberger-Janssen conjecture is correct, we find 
strong indications that 82 ~ 28, by using the corrected density. There is a small 
dependence of 8 on the diffusion factor d. We attribute this to the differences in 
the amplitudes of the second- and higher-order corrections. This induces a (small) 
systematic error in the values for 8 given in Table 1. 
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